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Abstract 

We propose a method, Temperature Integration, which allows an efficient calculation of free energy differences between 
two systems of interest, with the same degrees of freedom, which may have rough energy landscapes. The method is based 
on calculating, for each single system, the difference between the values of InZ at two temperatures, using a Parallel 
Tempering procedure. If our two systems of interest have the same phase space volume, they have the same values of In Z 
at high-T, and we can obtain the free energy difference between them, using the two single-system calculations described 
above. If the phase space volume of a system is known, our method can be used to calculate its absolute (versus relative) 
free energy as well. We apply our method and demonstrate its efficiency on a "toy model" of hard rods on a 1-dimensional 
ring. 



1 Introduction 

n 

Calculating free energy differences between two physi- 
cal systems, or between two thermodynamic states of the 
same system, is a topic of considerable current interest. 
The problem arises mainly in soft condensed matter, es- 
pecially in studies of macromolecules such as proteins or 
RNA. When the systems in question have complex en- 
ergy landscapes with many local minima, generating an 
equilibrium ensemble of configurations in reasonable run- 
ning time becomes a major challenge for computational 
physics. Indeed, a variety of advanced methods and al- 
gorithms have been introduced to answer the challenge, 
both in the context of Molecular Dynamics and Monte 
Carlo (for recent reviews see (TJ [2J [3] ) . 

Free energy difference between two systems can be cal- 
culated using equilibrium methods (as used by us) and 
non equilibrium methods. Existing equilibrium methods 
are composed of 3 stages: (1) selection of intermediates 
that interpolate between the systems (2) ergodic sampling 
of the system at each intermediate and (3) calculation of 
the free energy difference between the systems using one 
of the methods mentioned below. The commonly used 
methods include Bennett Acceptance Ratio [4], Weighted 
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Histogram Analysis Method [3] Exponential Averaging / 
Free Energy Perturbation [6] and Thermodynamic Inte- 
gration (Thi) HCOIB]. 

Non equilibrium methods measure the work needed in 
the process of switching between the two Hamiltonians. 
These methods use Jarzynski relations [9] (fast growth is 
one of their variants |10j ) and its subsequent generaliza- 
tion by Crooks 

Free energy differences are calculated in several con- 
texts, including binding free energies [12] [13], free ener- 
gies of hydration [2], free energies of solvation [TS] and 
of transfer of a molecule from gas to solvent [3] . Binding 
free energy calculations are of high importance since they 
can be used for molecular docking |16j and have potential 
to play a role in drug discovery |17j . 

Most of the applications mentioned above can be tack- 
led from a different direction using methods which mea- 
sure the free energy as a function of a reaction coordinate. 
These methods include Adaptive Biasing Force [18] and 
Potential Mean Force [7]. 

Our novel method, Temperature Integration (Tempi) 
can, in principle, be used instead of equilibrium meth- 
ods. In order to demonstrate the advantages of Tempi, 
we chose to introduce the idea in the context of Ther- 
modynamic Integration (Thi) [5J [7J [5]. Thi is based on 
simulating a set of systems defined by different values of 
a parameter < A < 1, where the two systems that we 
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wish to compare are realized when A = or 1. The free 
energy difference is given as an integral over A, which is 
evaluated numerically. Hence L, the number of values of 
A one needs, depends on how fast the integrand varies, 
which in turn is determined by the dissimilarity of the 
two compared systems. In general, the optimal choice of 
the intermediate systems is a challenge |17) . 

Since in many cases of interest each of the systems stud- 
ied has a complex energy landscape with minima sepa- 
rated by large barriers, equilibration times arc long. A 
favored choice to alleviate this problem is Parallel Tem- 
pering (PT) e.g. [19] a variant of the replica exchange 
method [20] [21] . This technique necessitates equilibration 
of a system of N particles at a set of n ~ \/~N inverse 
temperatures /3fc,fc = 1, ...n (where N is the number of 
particles). 

Combination of ThI and PT has been suggested by oth- 
ers [22] [23] as an efficient way [24] to calculate free energies 
of such systems. Since simulations of n replicas of the sys- 
tem are performed at each of L values of A, using PT with 
ThI calls for simulations at a set of L x n points in the 
A,T plane (see Figure [TJ. 

Our novel method, Tempi, uses the temperature dimen- 
sion, explored by parallel tempering, for the calculation of 
free energy differences; In effect the replicas, simulated in 
the parallel tempering procedure, are used as intermedi- 
ates for the calculation of free energy differences. Thus, 
the need for sampling both T and A dimensions is elim- 
inated. Furthermore, since in Tempi the internal energy 
(H) is a monotonic function of /3, the choice of intermedi- 
ates is no longer a challenging problem [ 17] (see Appendix 
for details), and the calculation is much easier to verify. 

Temperature Integration is based on calculating, for 
each system, the difference between In Z at the tempera- 
ture of interest and at a high temperature, using parallel 
tempering procedures. In case the two compared systems 
have the same phase space and hence Z at high-T, the dif- 
ference between In Z of the systems at the temperature of 
interest can be calculated and the free energy difference is 
obtained. For cases when the two compared systems have 
different Z value at high-T, we use an additional method- 
ological advance that enables the comparison. For some 
systems the calculation of absolute free energy values, us- 
ing Temperature Integration, is also feasible. 

The structure of the paper is as follows. In Section [231 
we introduce the method of Temperature Integration. In 
Section [3] wc describe how to compare two systems with 
different values of partition functions at high tempera- 
tures. In Scction|4]wc apply the method to a toy problem, 
demonstrate a calculation of absolute free energy values, 
and compare its performance with that of ThI combined 
with PT. The work is summarized in Section [5] 



2 Calculation of AF by Tempera- 
ture Integration 

Consider two systems, denoted by A and B, at a given 
temperature T\ = , between which we want to calcu- 
late the free energy difference AFa->b {Pi)- We will first 
assume that the two systems have the same degrees of 
freedom and phase space volume - so as f3 — > they have 
the same value of the partition function (the assumption 
of having the same value of the partition function will be 
relaxed in section [3]) . One of the most commonly used 
methods for calculating such free energy difference be- 
tween two such systems is Thermodynamic Integration 
[2J[7J[5], which we now briefly describe. 

2.1 Thermodynamic Integration 

Denote the Hamiltonians of the two systems by H A (Ct) 
and H B (n), where ft denotes the coordinates of the sys- 
tem. Noting that the two systems have the same coordi- 
nate space, we define a A-weighted hybrid system, charac- 
terized by the Hamiltonian H(\): 

h(x, n) = \H B (n) + (i - \)H A (n) (i) 

As shown in [3] [5] , the free energy difference is given by 

l 

AF A ^ B (M = J /^f\dA= (2) 
o 

i 

J[(H B (n)) x -{H A (il)) x ]dX (3) 
o 

where (X) x denotes the equilibrium average of X in 
the ensemble characterized by H(X). The expression for 
AFa-+b (A-)i written explicitly, takes the form: 

AF a ^bWi)= (4) 

} r [H B (n)-H A {a)]e-^ XHa ^+o--^mdn,. 

J J W) A 



(5) 

with 

Z{X) = J e -WAffB(fi)+(i-A)ffA(n)] dn ( 6 ) 

The integration is performed numerically, with the inte- 
grand evaluated at each one of a set of values of A by 
Monte Carlo simulations. As implied by ([5]), the two sys- 
tems are required to have the same degrees of freedom. 
The complexity of ThI is proportional to the number of A 
values I/, required to estimate the integral within a given 
error range. In the best case scenario of (H B — H A ) X 
being a monotonic function of A, this number increases 
roughly linearly with AF. Hence L is roughly propor- 
tional to the number of particles: L ~ AF ~ N. 
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2.2 Thermodynamic Integration with 
Parallel Tempering 

In many cases of interest both systems A and B, and hence 
also all the A- weighted intermediate systems, have rugged 
energy landscapes with many local minima. As the decor- 
relation times grow exponentially with AE/k^T , where 
AE is the energy barrier between nearby valleys, equili- 
bration times in these systems can be rather long. In or- 
der to overcome this problem, and obtain the equilibrated 
thermodynamic averages (H(Q)) X , one can use the Paral- 
lel Tempering procedure |25l 126] . However, implementing 
both thermodynamic integration and parallel tempering 
yields an unnecessary overhead in running time, as we 
will demonstrate. 

In order to implement thermodynamic integration we 
first choose a set of values Xi, i = 1, ...L, that will enable 
us to have a good sampling of the function (H(tt)) x for 
the integration in equation ([3]). Thus, L is related to the 
desired precision of the integration. 

In principle parallel tempering should be performed for 
each Aj, simulating each of the m A^-weighted systems 
over a set of temperatures, given by Pk, k = 1, ...n ~ y/~N 
[25] . Finally, using the calculated values of (H(fl)) x . at 
a temperature of interest Pi , we approximate the integral 
of Equation <j3j> by a sum over L terms, to get the free 
energy difference AFa-*b (Pi)- 

The procedure involves running Monte-Carlo simula- 
tions over L A-values, and n /3-values, that is, over a grid 
oi L x n instances of the hybrid system. An illustration 
of this grid is presented in figure [TJ 
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Figure 1 : Illustration of the grid of values over which the 
Monte-Carlo simulations are performed in the two proce- 
dures: ThI with PT and Tempi. The lowest values of p 
correspond to some very high finite temperature. 



2.3 Calculation of AF by Temperature 
Integration 

We present now our method, which obtains the free en- 
ergy difference AFa^b (Pi), using only simulations that 
are done in the process of parallel tempering, performed 
for the two systems A and B, (eliminating the need for 
simulations at a set of Xi values). This method can be 
applied to any two systems that have the same degrees of 
freedom f2. 

As P — > 0, the limiting value of the partition function of 
a system yields the phase space volume. In particular, if 
systems A and B, which have the same coordinate space 
{i~2}, have the same j3 — > limit, we have 



Z B (0 -> 0) = Z A (p 



o) = J <m . 



(T) 



In this case we can use the following identity to obtain, 
for a finite pi, the difference of the free energies: 



\nZ(pi)-\nZ(p^Q) = ^-dp = - (H) dp 

Jo d P Jo 

(8) 

Using equations (J7J and © we obtain 



AFa^b(Pi) = 7T I ln ^4 (Pi)- In Z B {fix)] 
Pi 



f 

Ji 



01 Pi 

(H B ) dp- J (H A ) dp 
o 



.(9) 



For each of the systems A and B we estimate the integrals 
on the right hand side by procedures of parallel tempering, 
sampling the system at a series of values pi, . . . , p n . We 
choose values such that the highest temperature sampled 
(corresponding to /3„) is much larger than the internal 
energy of the system at p\, so Z(P n ) ~ Z(p — > 0) is 
satisfied. 

3 Comparing two systems with 
different partition functions at 
high Temperatures by introduc- 
ing a cutoff 

The condition of Equation ([7J poses a problem for any two 
systems that have different partition functions values at 
high -T. In order to satisfy the condition stated above we 
had to set a cutoff over the interactions, E cuto s ■ We show 
that our results do not depend on the choice of E cuto s 
and P n , as long as -©cutoff is much larger than any typical 
interaction energy in the system at Pi , and /3~ 1 ^ -E cu toff ■ 
Note that this use of cutoff is general and can be used for 
any energy term that differentiates between the systems 
at the high T limit. 
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The proposed calculation of the free energy difference 
between the two systems at the temperature of interest 
Px is legitimate only if our choice of the cutoff energy 
has a negligible effect on the partition function value of 
each of the two systems at /3i. In addition, the highest 
temperature used, corresponding to /3„, must be such that 
the equality of the partition functions of the two systems 
is satisfied to a good accuracy. 

We denote the Hamiltonian with the cutoff energy by 
H' , and write the requirements stated above explicitly as 
follows: 

lnZ B 09i , H) ~ lnZ B (ft , H') (10) 

lnZ A (f3x,H)^lnZ A (f3x,H') (11) 

lnZ B (p n ,H')~lnZ A (p n ,H') (12) 

In order for the cutoff to have a negligible effect on the 
partition functions at the temperature of interest it has to 
be set to a value that satisfies 



E, 



cutoff 



> kaT x . 



As for j3 n , if the cutoff energy satisfies 
-Ecutoff -C ksT n , 



(13) 



(14) 



the systems will have almost equal probability to be in all 
the regions of their phase space, including ones which were 
restricted due to high energy values. Thus, the partition 
functions values of the two systems will be almost equal. 
Hence if these requirements are satisfied one can write: 

lnZ A (p u H)-lnZ B (Pi,H) ~ 
lnZ A (f3x,H')-lnZ B ([3x,H')~ 

lnZ B (p n ,H')-lnZ B (px,H') 
-[lnZ A {p n ,H')-lnZ A {px,H')] (15) 

Using the identity in Eq. ([8]), we can write: 

AF A ^ B (ft ,H) = fc [lnZ A ( ft, H) - lnZ B (ft , H)] ~ 



0i 



J^(H' B )d$-J"UH' A )df3 



(16) 



So the calculation of free energy will be negligibly affected 
by the use of a cutoff energy as long as we fulfill the rele- 
vant conditions. 

This use of cutoff energy is relevant also to ThI and 
similar methods. Consider for example two systems, Hq 
and Hx that have different steric constraints (resulting 
in different j3 — > limits), Then e.g. at A = there 
may be micro-states with finite statistical weight e~^ H " 
and infinite energy in Hi , so the sampling of the internal 
energy (Hx)x=o is infeasible. 

It can be seen that since the partition function of a 
system that has the Hamiltonian with the cutoff H' is 
almost equal to the one of a system with the Hamiltonian 
H at fix. Thus 



AF A ^ B (ft, if) ~ -j- [lnZ A (j3x,H') 
Pi 



lnZ B (/3x,H')}, 
(17) 



and ThI can be implemented for systems with H' and 
yield almost the same result. 

In conclusion, with the use of cutoff energy in ThI, that 
enables us to sample the integrand in all cases, the free 
energy difference between any two systems that have the 
same degrees of freedom can be calculated. 

4 Applying Temperature Integra- 
tion to a toy model 

In the following sections we present the results obtained by 
applying Temperature Integration with interaction cutoff 
to a "toy model" , of N = 8 — 25 particles on a ring (one 
dimension, periodic boundary conditions) in an external 
potential, interacting via a hard core potential. For this 
model we evaluated the free energy difference between two 
systems, which differ in the size of their hard cores. We 
did this in two ways: First we performed ThI at a set of A 
values, where at each value the corresponding system was 
equilibrated using PT; second, we used Temperature Inte- 
gration. We compared the results as well as the amount of 
(computational) work needed, for the two ways, to achieve 
similar accuracy. We estimated the gain in work (number 
of Monte Carlo steps needed) and the way it scales with 
the number of particles. 

4.1 Definition of the model 

Here we demonstrate the method for a toy model in one 
dimension with periodic boundary conditions: we place N 
particles on the unit circle, with the position of a particle 
defined by an angle (9 = [0, 2tt]). The particles are in an 
external potential, given by 



V (0i) = V ■ cos (26,) 
and have "hard core" interaction: 



(18) 



TTfO o \ - / 00 dist(6i,0j) <W nQ , 
U V*> e i)-\ diatpM^W (19) 

Here W is the size of a particle's "core", and dist (9i,9j) 
is the angle difference between the particles. 

In order to apply Tel we introduced a cutoff of this 
potential, replacing U (9.^9j) by: 



U 



i) 



UcutoS dist (9i 
dist (9, 



9 3 ) < W 
9j) ^ W 



(20) 



where t/cutoff is the energetic cost of two particles having 
overlapping cores. The effect of the cutoff on the results 
of the calculation is negligible as shown below in 14.51 

We calculate the free energy difference between two sys- 
tems A and B, that have the same number of particles and 
different values of W. Specifically, we set the constants 
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Vo, W in the systems A, B to have the following values]: 



V, 



OA 



V, 



OB 



w A 



7T 1 

4 TV' 



We 



7T 1 

8 N 



(21) 



We work at temperature fee?! = 0.6157 so that the 
barrier associated with the external potential is signifi- 
cantly higher than k B Tx, and set f7 cu toff = 70, so that 
^cutoff ^ k B Ti and hence condition (flU]) is satisfied. 



4.2 Details of the Monte Carlo Simula- 
tion 

In the initial configuration all the particles were placed in 
the interval [0,7r], at equal distances. 

The local Monte Carlo move consists of randomly choos- 
ing one of the particles and changing its position 6 to 
6 + A6, where A8 £ [~ ifg, ifg] is selected with uniform 
probability. 

All the Monte Carlo simulations are performed using 
Parallel Tempering; that is, after a certain number of local 
moves (500 in our case), we attempt to exchange config- 
urations between systems at adjacent temperatures. We 
simulated systems with 8,10,12,14,17,20 and 25 particles 
and performed up to (1, 1.1, 1.2, 1.35, 1.5,2)-10 8 MC moves 
in total for each one respectively. 

4.3 Using Thermodynamic Integration 
and PT 

In order to compare the results obtained by Tel we used 
ThI in combination with PT. The set of temperatures at 
which we worked was selected as follows. 

The highest temperature was chosen as 3Vb/fcB, to en- 
able particles to cross easily the energy barrier (2Vo) of 
the external potential. 

The temperatures for both systems A and B were first 
selected so that the acceptance rates for exchanging con- 
figurations at neighboring temperatures will be between 
0.25 and 0.35. Then, the set of temperatures to be used 
was of that system in which the product of the accep- 
tance rates was lower (i.e using the set of temperatures 
with higher density). 

The integral in (fTI)]) was calculated by adding sampling 
points according to global adaptive Simpson's quadrature 
(GASQ) [271 [21] • That is, for each value of A, a simula- 
tion of a compound system was performed at the selected 
set of temperatures, and (Hb{^)) x r 1 — (Ha{^))\ p 1 was 
calculated. 

At each division of an interval, we sampled the internal 
energy at another 4 A values (according to GASQ [2"Tll2"8] ). 
calculated the integral with the current set of intervals and 
displayed the result as a function of the total number of 
MC steps performed. 



The -j* factor was added to maintain constant density of par- 



4.4 Using Temperature Integration 

For each system, we evaluated numerically the integral: 



fc B Ti / (H)d(3 



(22) 



We chose k B T n to be given by k B T„ = 400 *U *N*(N — 
l)/2 in order to satisfy condition (|T4")) . 

This was done iteratively. In each iteration we per- 
formed parallel tempering on a different set of tempera- 
tures and calculated the internal energies. Then, we added 
the results to those of the former iteration and calculated 
the integral (f2"2"|) numerically according to the current set 
of sampling points (j3- values and the corresponding inter- 
nal energies). 

After each calculation of the integral, we registered the 
result as well as the number of Monte Carlo steps per- 
formed in total until a stop criterion was reached (see 
Appendix). 

The method by which we chose the sets of /3- values was 
based on the global adaptive Simpson's quadrature and 
was suited to maintain optimal acceptance rates between 
the systems (as explained in the Appendix). 

4.5 Absolute Values of free energy and 
Verification of the method 

The highest temperature T n was chosen to be much higher 
than the total interaction energy, so the partition func- 
tion becomes the phase space volume, tt. Hence, when 
the phase space volume of a system is known at T n , we 
can calculate the absolute (not relative) free energy of the 
system: 



F = -k B T x 



n 



ln 73lv+ / WW 



Pn 



(23) 



where the integral in equation (|23[) can be estimated using 
the methods described in Section [231 

In many systems (including our toy model) when we 
neglect all interactions (including steric) the phase space 
volume is known and the values of free energy can be 
calculated according to (|2"5)l . enabling the comparison of 
free energies of systems with different degrees of freedom. 

To assess the accuracy of Tempi (also in the context of 
absolute free energy values) for our toy model, we com- 
pared its result (for N = 3 particles) with "exact enumer- 
ation" numerical evaluation of the free energy, obtained 
by performing with high precision the integral 



F = -feeTi In 



fe-^ H dn 



(24) 



tides. 



In this integration, the microstates in which the steric lim- 
itations were violated weren't taken into account, enabling 
us a liable comparison of the results according to 
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4.6 Results 



ThI (with PT) and Tempi converged to the same asymp- 
totic value of the AFa^b after a large enough number MC 
steps. The number of MC steps required for each method 
to approach the asymptotic value is, however, very dif- 
ferent. In Fig. [3 we present the relation between the 
values of AFa^-b, calculated by the two methods, and its 
asymptotic value, as a function of the total number of 
the MC steps used in the calculation, for systems with 8, 
17 and 25 particles. Evidently, in the Tempi procedure 
the convergence of AFa^b to the asymptotic value was 
significantly faster than in ThI. 





ThI 

Tempi 










6 ' 



1.1 ■ 



N = 17 



V 



N=25 



x10 10 x10 10 

Monte Carlo Steps (MCS) 

Figure 2: The result for AF, normalized by the "asymp- 
totic" value, as a function of the number MC steps using 
ThI and Tempi, for systems with 8,17,25 particles. 

Another measure of the relative efficiencies of the two 
methods is the number of MC steps needed to achieve a 
desired level of accuracy. In Fig. [3]we present the number 
of MC steps needed to approach the asymptotic value to 
within 1%, by both methods, as a function of the number 
of particles. The ratio between the number of MC steps 
needed for convergence in ThI and in Tempi increases 
steeply with particle number, as demonstrated in Fig. [3] 
We also computed the absolute free energies for N = 3 
particles, by both Tempi and by exact enumeration. The 
numerical results obtained by Tempi, using (|23|) . were: 



Fa 



Inn (T„) + Jf (H) dj3 



= -k B Ti [5.51 + 12.72] = -11.23 
F B = -k B T! [5.51 + 13.40] = -11.64 

where the phase space volume is = (2ir) . 



(25) 
(26) 
(27) 

Numerical 
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Figure 3: MC steps needed for 1% convergence in both 
methods as a function of the number of particles, N. 



integration over the whole phase space yielded the same 
values, confirming the validity of our method in general 
and in the context of absolute free energy values. 



5 Discussion 

We presented Temperature Integration, a method to cal- 
culate the free energy difference between two systems at 
some inverse temperature f3\ . Temperature Integration is 
an efficient method since the temperatures used in the 
parallel tempering procedure are used as intermediates 
in the calculation of free energy. Moreover, the method 
for choosing the intermediates for the calculation of free 
energy in Tempi is general (see Appendix for details). 
Hence, the method is robust, which is important for au- 
tomation and high-throughput use. 

In the calculation, we performed parallel tempering pro- 
cedures for each of the two systems over sets of temper- 
atures between /3i and /3„, where /3" 1 -Ecutoff - the 
maximal interaction energy of the system. Then, we use 
the internal energy values obtained at each temperature 
to estimate numerically the two integrals in equation ([9]), 
and hence the free energy difference between the two sys- 
tems. 

Furthermore, absolute values of the free energy can be 
calculated for systems in which the phase space volume is 
known (when all interactions are neglected) as stated in 
section 14.51 

Tempi can be used also for systems with a smooth en- 
ergy landscape (where PT is not needed). Tempi has the 
advantage of simplicity (saves programming time) since 
the simulations are performed only on the two original 
systems (i.e. using only the two "pure" Hamiltonians). 

While in other methods for calculating free energy dif- 
ference between two systems the choice of appropriate in- 
termediates remains a challenge |17| , in Tempi the internal 
energy (H) is a monotonic function of /? and as a result 
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the intermediates can be easily chosen (see appendix for 
details) and the calculation is much easier to verify. More- 
over, the monotonicity of the function may result in less 
intermediates and facilitate the calculation (the number 
of integration points scales as the free energy difference 
which is roughly linear with the number of particles N). 

The method was applied to calculate free energy differ- 
ences in a toy model of hard rods on a 1-dimensional ring. 
It was also applied to calculate free energy differences be- 
tween interior loops of RNA that differ in a base (See [29] 
for details). 

There is a provisional patent pending that includes the 
contents of this paper. This work was partially supported 
by the Leir Charitable Foundation (ED, AF) the Einstein 
Center for Theoretical Physics (NC) and the National Sci- 
ence Foundation under CHE-0713981 (CHM). We want to 
acknowledge Amir Marcovitz for his assistance. 

Appendix: Integration method in 
the Temperature Integration proce- 
dure 

We introduce the method in which the sets of /3-values 
are chosen in the process of sampling the function in the 
Tempi procedure. These sets have to be chosen in a way 
that will minimize the total error calculated according to 
the Simpson's method and will satisfy optimal acceptance 
rates in the PT procedures. 

We chose a temperature Tpx - higher than the energy 
barriers in the system. In each PT procedure we kept 
constant the number of /3-values sampled within the range 
[Ti,Tpt], used in the PT procedure, and added them to 
a set of /3-values in the range [Tpt, T„]. 

In this method we first sample over the set of /3-values 
chosen to satisfy optimal acceptance rates in the system, 
including the points (/3„ + /3pt) /2 and j3 n , as defined here 
and in Sec. l4~4l 

Second, we sample over a set of /3-values that bisect the 
intervals of the former set. In the subsequent bisecting we 
take into account the facts that in the range [Ti , Tpx] the 
requirements for PT go together with the ones needed for 
optimal sampling of the integral, and that in the range 
[TpT,T n ] the two requirements don't necessarily correlate 
but the temperatures can be sampled independently of the 
value of the other temperatures in the set. 

The 5 sampling points in the range [TpT,T n ] we now 
have form a subinterval according to ASQ. We further 
sample this subinterval according to GASQ until the max- 
imal error in the range (defined by AxjAj/i) is smaller 
than the maximal error in the range [Ti , Tpt] ■ The group 
of subintervals in this range will form a vector which will 
be called u h i g h temps- 

We now further bisect the sampling points in the range 
[Ti,Tpx] in 2 PT procedures, with the subinterval with 
the maximal error in «hi g h temps (according to GASQ). We 



repeat this step, necessitating this time 4 PT procedures 
and the sampling of the 2 subintervals with maximal error 

m Vhigh temps- 

Then, we generate 2 arrays of vectors of subintervals in 
the first range. The first consists of the odd subintervals 
and the second of the even subintervals, each vector con- 
sisting of the subintervals that will be derived from the 
original subinterval. We perform a loop in which in each 
iteration, we choose one of the 2 arrays. Then we select in 
each vector in the array the subinterval with the maximal 
error and we add to the chosen set the subinterval with 
the maximal error in Whigh temps- We bisect the chosen set 
of subintervals and since each bisecting includes 4 sam- 
pling points, we perform 4 PT procedures in order to do 
so. Then the bisected subintervals are placed instead of 
the subinterval from which they were derived, and the in- 
tegral is calculated. Thus, we proceed until a certain total 
error or to a maximal number of iterations is reached. In 
each calculation of the integral, the number of MC steps 
performed is registered. 
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